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1. Introduction 


The primary model for analyzing the dynamics of many species' populations is the predator-prey model, 
which was developed by Lotka and Volterra [1, 2]. It has a broad range of applications in many different 
scientific fields, including chemical processes [3, 4], bioparticle granulation [5], and the interaction of 
microorganisms and ecosystems [6, 7]. The Lotka-Volterra type predator-prey system has been the focus of 
a lot of research in recent years. Zhu and Yin focused on the competitive Lotka-Volterra model in random 
situations [8]. Badri et al. focus on installing Lotka-Volterra systems, a particular family of nonlinear 
quadratic systems, and possible balance points [9]. This paper's objective is to numerically study the 
predator-prey models controlled by the following system of nonlinear ordinary differential equations 


= = f (x(t), (0), a1 (0), -.-, an) 
& = g(x(t), y(t), br(t), » bat) 


where the positive functions a;(t), b; (t) typically give relative measurements of the effect of dimensional 
parameters [10] and x(t), y(t), respectively, reflect the population densities of prey and predator at time t. 
The numerical solution of system (1) is not straightforward due to the structure of the functions f and g, 
hence it is essential to create trustworthy numerical algorithms. Numerous articles have addressed the 


() 
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numerical solution of predator-prey models to compare the effectiveness and reliability of various 
approaches numerically. For instance, the Adomian technique has been quantitatively tested using a 
predator-prey model [11-14]. Additionally, He's variational approach was researched and used to simulate a 
predator-prey relationship [15]. The predator-prey model has also been used using nonstandard finite 
difference techniques [16]. The differential transformation method (DTM) was used to analyze a predator- 
prey model with constant coefficients over a brief time horizon [17]. The ratio-dependent predator-prey 
system with continuous effort harvesting was approximated using the homotopy perturbation approach [18]. 
The nonlinear imprecise prey-predator model with stability analysis was implemented using the homotopy 
perturbation and variational iteration methods [19]. Many researchers have presented studies on the 
characteristics of prey and predator systems such as coexistence, stability and extinction [20-22]. A 
numerical study is presented in a ratio-dependent predatory system disturbed by time noise. The base model 
was calculated using two charts. Euler forward diagram and finite difference diagram [23]. Using a Holling 
type-II predator functional response, a reaction-diffusion-advection predator-prey model was shown to be 
stable [24]. A multispecies ecological and epidemiological mathematical model was designed for scenarios 
where interacting species compete for the same food sources and the prey are infected. The stability of the 
aforementioned method is examined by the Von Neumann stability analysis [25]. What was presented above 
reflects the importance of studying the system and the importance of treating it with different simulation 
methods. In general, in many cases, it is difficult to find a solution to nonlinear differential equations using 
the methods of integrative transformations due to the nonlinear terms in them. Moreover, according to our 
information, the process of merging analytical methods and integrative transformations leads to a reduction 
in the number of computational operations and reduces the difficulty of analytical methods are used alone. 
Many studies have focused on this area [26-29]. Many analytical methods used to solve differential 
equations, and one of these methods is Akbari-Ganji’s method developed by Ganji and Akbari and could be 
used to solve a variety of nonlinear ordinary and partial differential equations. It is characterized by its ease 
of application to find solutions [30]. In addition, Pade’s approximation is the best approximation of a 
function near a certain point by a logical function of a specific order. It improves accuracy and the 
convergence of solutions by improving the field, which was developed by Henri Pade' in 1890 [31]. A recent 
integrative transformation is the Shehu transform, which was developed by Maitama and Zhao in 2019 [32]. 
And the fact that these aforementioned methods have not been used before to solve the prey and predator 
systems. Consequently, this has prompted us to propose a new method. This new method combines the 
algorithms of the Shehu transform method, Akbari- Ganji’s method, and Pade’s approximation method 
(SAGPM). This paper involves seven sections, of which this introduction is the first. In Section 2, We 
present the proposed technical algorithm, and we applied it to different systems of predator- prey models 
(Section 3). Next, we introduce the numerical results and compare these results with other works (Section 4). 
In Sections 5 and 6, we study the convergence and stability analysis of predator-prey systems. The last 
Section summarizes the major findings of this study. When applying the proposed approach and comparing 
it with the previous works that dealt with the same problem, the results have shown the effectiveness and 
efficiency of the proposed approach in terms of high accuracy and good convergence. 


2. The new SAGPM Algorithm 


The basic idea of the SAGPM is based on the Shehu transform method and Akbari-Ganji’s method with 
Padé approximate algorithms, which will be mentioned in this section. 


2.1. Shehu transformation method [32]: 


Shehu transform, a new integral transformation that Shehu and Zhao introduced in 2019, is a generalization 
of the Laplace and Sumudu integral transforms. The authors have used it to resolve both the ordinary and 
partial differential equations. 


The Shehu transform is found across set A: 


A={f(t):43 N,n4,nz > 0,|f()| < N exp (), if te (-1)/ x [0, 00) } 
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by the following integral 

Sift)] = F(v,w) = J” exp) fit) at 

=lima-soo J, exp(—) F(t) dt, v>0,u>0. (2) 
The inverse Shehu transform is given by 

S/F(v,u)] =f(t), fort> 0. 

Equivalently 

ft)= SUF (vw) for + exp) Fru) dv, (3) 


where v and wu are the variables of the Shehu transform, and @ is the real constant and the integral in Eq. (3) is 
taken along v = a in the complex plane v = x + iy. 


2.1.1. Shehu transformation of derivatives and some functions 


If S{f(t)} = F(v, wu) then 
1) S{4} =? Fu) - FO 


2) s{f(} = P@w — yaa (2)? pw 
3) Si}=5 

4) S{exp(at)} = — 

5) S{cos(t)} = 


vu 


6) S{sin(t)} = ae 


7) s{-} = er n= 0,1,2,... 


nN 
2.2. Akbari-Ganji’s method (AGM) 


This method, which was developed by Akbari and Ganji, is an excellent computation methodology that may 
be utilized to solve various nonlinear differential equations. In this method, the solution is taken to be a finite 
series, hence the answer is obtained by resolving a sequence of algebraic problems. 


The differential equation for a function x(t) and its derivatives can be written as follows in order to use 
AGM: 


In equation p, = f (x, x',x",...,x%) =0, x =x (t) (4) 

Dx is the nonlinear differential equation of m“"order derivatives with boundary conditions: 

At t=0, x(t)=X, x (t= X1,..., XD) = X17 (5a) 
AKL, 4) =e 5% WH tient” O=%,., (5b) 


The differential equation's solution is taken into consideration to solve Eq. (4) in relation to the conditions 
(5a and 5b); 


x(t)=Lixo a tt (6) 


Eq. (6) can be solved with high accuracy by selecting more terms. If the series (6) has (n) degrees, then there 
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are (n + 1) unknown coefficients that must be determined to find the solution to the differential Eq. (4). For 
Eq. (6), the boundary conditions (5a) and (5b) are applied as follows: 


We have the following at t=0. 
x(0) = ay = Xp 


x'(0) =a, =x, 
7 
x" (0) =a, =x, M 


when t=L 


x(L)=a,+a,L+...+a,L' =x, 
x'(L)=a,+2a, L+...+na, =x, (8) 


x"(L) = 2a, + 6a, L...tn(n-la, LD? =x, 


After substituting Eq. (6) into Eq. (4), and after applying the boundary conditions to it, we obtain: 


Dy = f(x(0),x'(0),x""(0),....x°" (0) 
Dy = f (XCD), x"(L), "(Dyes 8" (L)) (9) 


Application of the boundary conditions on the derivatives of the differential Eq. (9) is: 


sie ellie gl tis 
DFE) x (LE) 21D) gees 2 (D)) 


" 


Lok. eee (11) 
"FC BP aya) 


Equations from (7) to (11), (n + 1) equations may be worked out, and so (n + 1) unknown coefficients of 
Eq. (6), such as ay, Q,, Qz,..., and a, can be computed. By locating the coefficients of Eq. (6), the solution 
to Eq. (3) will be achieved. 


2.3. Padé approximant 


George Frobenius, who introduced the concept and researched the characteristics of the rational power 
approximation, is credited with creating the Padé approximant, a specific and classical sort of rational 
approximation. By expressing it as the quotient of two polynomials with various degrees, Henri Padé made 
substantial contributions about 1890. In comparison to truncating the Taylor series, it provides a better 
approximation of the function, especially in cases where the Taylor series does not converge and contains 
poles. This makes it superior to the Taylor series expansion. The approximation has been widely used in 
computer science to determine time delays[33-36]. 


Padé approximation is a ratio of two polynomials that come from the Tayler series expansion of a function 
x(t) and define as [31] 


L n 
Pp _ Yn=0 Ant (12) 


> Se baat 
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where bh = 1 . 


The function x(t) written by 


x(t) — ar G.t° a 
thus, 
yop ant” 
nao Cnt = on (14) 
or 
ee ee oe Oe 
' : . 1+bytt+bot?+- 


from Eq. (14), obtain the following system equations 
Ag = Co 
ay = Cy + Cob 


az = C2 + C,b,; + Cob 


solve the above system of equations for a, and b, let the numerator degree to | and the denominator 
degree to m. 


Ao = Co 
ay = Cy + Cob; 


Az = Co + C1), + Cobz 


a, = C, + C154 + Ci-2b2 fee Cob) 
0 = C44 + (by + Cy bz + ++ C}_-m4ibm 


O = Cig + Cr41d1 + Cybz +++ + Ci-m+2Dm 


O = Clam + Ci4m—191 + Ci4m—2b2 + + + Dm 


Therefore; the fundamental notion of the new technique SAGPM for system (1) with initial conditions x(0) 
=a ,y(0) = f is summarized according to the following steps: 
Step 1: Taking Shehu transformation on both sides of (1), to get: 
d 
SCZ) =SEO,y(),a1(0,.-an(0)), 
d 
S()=S(g(x(0.y(),b1(0),---ba())), 


using the differentiation property of Shehu transformation and the above initial conditions, we have: 


X(u,v) = ——+—— S(f(x,y(0),a1(),..an(©)), 


Yu, v) ==" += S(g(x(,y(0).b1),..-bn(t))), (15) 


Vv 
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Step 2: Applying the inverse Shehu transformation on both sides of (15), to find out: 
x(t) =a + S7(—S(F(x(0,y(,a1(1),-.-an())), 


yt) = B + S*(—S(g(x(t),y(t), by, bn (t)))), (16) 
Step 3: Considering by the AGM as polynomials series with constant coefficients as follows: 
x(t)= Dikoai th » yt) = Lio bi t! (17) 


After substituting (17) into (3), Eq.(16) becomes as follows: 
dito ai t=a+ SSL Eo aj ie a bj thay san) 


Tho bi th =B + S*(—S(g(Dho a t!, Lik bi t! bi ),...Pn(2)))), (18) 
We apply the boundary conditions to obtain some values of the coefficients, and by continuing the 


derivation of the Eq. (18) and substituting the boundary conditions, we get the rest of the values a; and b;. 


Step 4: Applying the Padé approximation of an order [l /m] on a power series solution obtained by using 
(SAGPM). The values / and m are arbitrarily selected. In this stage, we get the final solution. 


3. Numerical experiments 


In this part, three different systems corresponding to predator-prey models are solved using the new 
technique. 


3.1. Example 1: 
The first model depicts the problem of some rabbits and foxes coexisting, with foxes preying on the bunnies 
and rabbits preying on clover, and with an increase and reduction in the number of foxes and rabbits, 


respectively [37]. The ordinary differential equation system shown below serves as the model's analytical 
representation: 


4 (€) = a x(t) — ay x1(1)x2(0) , x, (0) = 2 (19) 
Kp (t) = — By x2(0) + Bz x1 (1)x2(0) , x2(0) = 2 
with 
a, =a, = —t, Bi =f, = t. 
The exact solutions are given: 
14() = 2, m(th= —z 
2-e2 2-e2 
By taking Shehu transformation on both sides of (19), we get 
Kus) = SF + 2 S(ar a) - a2 1 Ox2(0) ; 
X,(u,s) = <e “lr = S(-B X2(t) + Bz x,(0)x2(t)) = 
Taking the inverse Shehu transformation on both side of (20), so, we get, 
X(t) = 2+S7*( - SC ay X(t) — a2 x1(t)x2(t) )) (21) 


X,(t) = 2+S7*( = S(— By X2(t) + Bz x1(t)x2(t) )) 
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By AGM , we must substitute (17) into (21), so, we get 
eo G th = 24+S71( - S( Oy Vika a th — a2 Vga; th * Vg bj t! ) 
Tobi th = 2+S*( = S(— Bi Lo bi t' + Bo Lio ai t' * Lo bi t')) 


When n =4, after simplification and offset values of a1, @2, B, and Bz, Eq. (22) becomes: 


(22) 


f (x10) = —t?agby + (—agb, — agb3)t® + (—agb4 — agb3 — agb,)t’ + (—ayb4 — azb3 — agb2 — 
ayb,)t® + (—b3a, — apby — agb, + (—by + 1)a4 — byag)t® + (—a,b2 — anb, + (—bp + Laz - 

b3a9)t* + (—a,b, + (—bp + 1) az — bzay + 40,)t? + ((—bo + 1)a, —apb, + 3a; )t? + (2a, + 

(—bp + 1)ay)t +a, =0 (23) 

g(x2(t)) = —t?ayb4 + (—azb, — agb3)t® + (—azbq — a3b3 — agb2)t” + (—ayb, — agb3 — agb2 — 
ayb,)t® + (—azb, — aby — b3a, + (—dg + 1)b, — boa,)t? + (—by az — a, bz + (—ag + 1)b3 - 

boa3)t* + (—a,by + (—ag + 1)b2 — bodg + 4b4)t? + ((—ap + 1), — ayby + 3b3)t? + (2b + 

(—dy + 1)bo)t + b, = 0 (24) 

The constant coefficients of Eqs. (23) and (24) which are {d,...,a, and bg,...,b4 } can be computed by 
applying boundary conditions in the form of: 


y(t =0)=2>a =2, x,(t=0)=2 >), =2. (25) 


On both Eqs. (23) and (24) and its derivatives, the initial conditions are applied as follows: 


f( we = 0)): a, = 0, Gee = 0)): b, =0 

fia) » ay by +49 + 2a, =0, 

g'(*2=9) : ay by + bo + 2b, = 0, 

fel10=O) ; — 2ay by— 2a, by +2af +aP= 6, 

g(1G=9). _ 24a, b, — 2a, by + 2b, + 6b; = 0, 

fia). — 6ay by — 6a, b, — 6a by + 6a, + 24a, = 0, 

g!"(1=%) . — 6a) by — 6a, by — 6ay by + 6b, + 24d, = 0. 
With the help of Maple software, unknown constant coefficients a; and b; can now be found by solving the 
above equations, these coefficients are illustrated below: 


a, = b, =0,a, = bz = 1,a3 = bz =Oanda, = by =3/4 
Then, x(t) =x,(t) = 2+ t? + 0.750000t* (26) 
Now, take the Padé approximant of Eq. (26) ,with | = 2, m= 2, we get the solutions of the system (19): 


2-0.5t2 
X41 12,2(©) = *2i2:31 (= 1—0.75t2 


3.2. Example 2: 


The second model takes into account the problem that the predator in the model is not significant from a 
commercial standpoint. The prey is continuously subjected to effort harvesting, which has no direct impact 
on the predator population. The availability of prey to the predator indirectly lowers the predator population. 
Furthermore, it is expected that the population of prey would rise simply logistically [38]. The following 
system represents this model: 


© 2023, CAJMTCS | CENTRAL ASIAN STUDIES www.centralasianstudies.org ISSN: 2660-5309 | 108 


CENTRAL ASIAN JOURNAL OF MATHEMATICAL THEORY AND COMPUTER SCIENCES _ Vol: 04 Issue: 08 | Aug 2023 


% (t) = x,(€)(1 — x, (¢)) — bz(t) - a 
X2 (t) = cz(t)- e x2(t) 
where 


(27) 


(t) = X(t) x2(t) 
x4 (t)+ x2(t)’ 


For the numerical simulations of the model (27), we take 


x,(0) = 0.5,x,(0) = 0.3,b = 0.8,c = 0.2,e = 0.5,r = 0.9. 


By taking Shehu transformation on both sides of (27), we get 


ee + * s(a(O(1- xO) — ba) — 10) 28 
84 4 ¥ Se 2(¢)- e x,(t)) ” 


X,(u,s)= 


X,(u,s)= 


Ss 


Taking the inverse Shehu transformation on both sides of (28), so, we get, 

x1 (t) = 0.5+S7'(— S(x,(O)(1 — x1(t)) — bz(t) — rx, (4))) i 

x,(t) = 0.3+ se Ge S(c z(t)- e x,(t))) 

By AGM, we must substitute (17) into (29), so, we get 

Vico a th = OS2kS*SCLihy a t'(1 — Dh a t= DZD = edo gist) 

30 

ob t' = 0.34 SS S( cZ(t) — e iL, b; t')) wa 

where 

At) = 


MES aj t! *yL9 bi t! 

He aj tty, bi ti- 
Whenn =3, after simplification and offset values of b,c,eandr, the Eq. (30) becomes: 
f(x.(0)) = Sta + 2taz + ay — (as + its + ta, + ay) (—t? a, va t? a> —, ta, a Ao + 1) + 


0.8 u + 0.9(t?a3 + t7a, + tay +a) = 0, (31) 


g(x2(t)) = 3t7b3 + 2th, + b, — 0.2 p+ 0.5(t3b3 + 2b, + tb, + by) =0, (32) 
Where 


_(t8a3+t?az+ta,+a9)(t3b3+t?b2+tb1+bo) 

~ £3(a3+b3)+t?(az+b2)+t(a,+b1)+apbo 
The constant coefficients of (31,32) which are {d,...,@3 and bg,...,b3 } can be computed by applying 
boundary conditions in the form of: 


x(t =0)=05->a, =0.5, x(t=0)=0.3 > bd) = 0.3 (33) 
On both Eqs. (31) and (32), and its derivatives, the initial conditions are applied as follows: 


0.8ap bo 
+0.9a, =0 
Ap + bo 


7( Milt = 0)): Qa, —Aj(1—ay) + 
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0.2agbo 
g( x(t = 0)): b, - aT 


fi(1G=0) : 2a, = a,(1 _ Ao) + Apa, + 0.8( 


abo i Agb, _ Apbo (ay + ba), 4 
Agtbo aot bo (ay + bo)? 
a, bo Agb, Ag bo (a, + a), +05b, =0 
. al — 


0.9a, =0 


Ag +by agt bo (dy + bo)? 
fi), 6a3 ae 205(1 = Ao) + 2a? + 2ap az 


g'(2=9) : 2b, — 0.2( 


406 azbo if a,b, _ aybo(a;b,) Agbz _ Apby(a, + by)  aAgbo(a, + b,)? 
“Vag tbo ap + do Ap + bo Ay + bo (dp + bo)? (ay + bo)? 
0.8apb9(2ag + 2b2) 


(a, Fb et 842 = 0 


gi(=0), 6b, 
= +( azbo a,b, aybo(a, + b,) Agbz Agby(a,+b,) agbo(a, + =~) 


Ay + bo 7 Ag + bo (ag + bo)” . Ag + bo (ao + bo)” (ap + bo)? 
Apbo (2a, + 2b,) 
0.2 ———————_ + 1.0b, = 0 
" (ag + bo)? . . 


With the help of Maple software, unknown constant coefficients a; and b; can now be found by solving 
equations, these coefficients are given below: 


a, = —0.35000, az = 0.19476562,a,; = —0.10728816, b; = —0.112500, b, = 0.01880859, and bz = 
—0.00112847 

Then 

x1(t) = 0.5 = —0.35000t + 0.19476562t? — ieee, Sa 3, 


34 
x, (t) = 0.3 — 0.112500t + 0.01880859¢? — 0.00112847¢3 eo 


Now, take the Padé approximant of Eq. (34) ,with 1 = 2, m= 1, we get the solutions of the system (27): 


_ 0.5 — 0.07457109077t + 0.001965388502t? 


th= 
X1 2,10) 1 + 0.5508578185t 

ne 0.3 — 0.09450058413t + 0.01205881279t2 
*2[2,41? = 1 + 0.05999805297t 


3.3. Example 3: 


The final model under consideration is the non-autonomous system of ordinary differential equations called 
a Lotka-Volterra model. In this model, time-varying values for the prey's growth rate, the predator's 
effectiveness in catching prey, the predator's mortality rate, and the predator's rate of growth are taken into 
account. It is significant to note that careful consideration must be given because this problem's coefficients 
are time-varying to acquire the proper recurrence equation system for the model. In other works [37, 39], the 
ordinary differential equation system shown below describes the aforementioned model: 


X(t) = a x(t) — a2 x4 (t)x2(t) 
Xp (t) =— By X2(0) + By x41 (t)x2(0) 


For the numerical simulations of the model (35), we take 


(35) 
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a, =4+ tan(t),az = exp(2t), B, = —2,B, = cos(t),x, (0) = —4and 
X2(0) = 4. The exact solution for these coefficients is 


x1(t) = — , XQ (t) = 4exp(—2t). 


By taking Shehu transformation on both sides of (35), we get 
4u u 
X,(u,s)= a a ets x1(t) — a2 x,(t)x2(t)) 
4 
X,(u,s) = = + = S(—By x(t) + Bo x1 (O)x2(0)) 


Taking the inverse Shehu transformation on both sides of (36), we get, 


(36) 


x(t) = —44+S$710¢ = SC a, x1 (t) — ay x4 (t)x2(t) ‘| 
. (37) 


X,(t) = 4+S7*( = S(— By X2(t) + Bz x1 (x2 (0) )) 
By AGM , we must substitute (17) in (37), so, we get 


io @ th = —4+S°1¢ = So Do a t' — a2 Vio a th * Viigo dj t! ‘ 
(38) 


Dito bi t! = 4+S*( = S(— Bi Lho bi t' + Bo Lko ai t' * Lh bi t' )) 
Whenn =2, after simplification and offset values of a@1,a@2,B, and B2, Eq. (38) becomes: 


f(x1(t)) = 2ta, + a, — (4+ tan(t))(t?a2 + ta, + ag) + e7*(t7a2 + ay + ag) (t?bz + thy + 
bo) = 0 (39) 


g(x2(t)) = 2tb, + by — 2t7b, — 2tb, — 2by — cos(t) (t?az + tay + ay) (t7bz + tb, + bo) = 0 (40) 


The constant coefficients of (39,40) which are {a ),Q,,Q2 and by,b;, bz } can be computed by applying 
boundary conditions in the form: 
y(t =0)=-4>-a, =-4, x(t=0)=4>bd, = 4. 
On both Eqs. (39,40) and its derivatives, the initial conditions are applied as follows: 

f( Kitt = 0)): Ao bo _ 4ao + ay = 0, 

g(x2(t = 0)): = Agbo — 2bo + b, = 0, 

fils) : 2ayby + agb, + a,b — ap — 4a, + 2a, = 0, 

g'\*2=)) : —ayb, — ayby — 2b, + 2b, = 0, 
With the help of Maple software, unknown constant coefficients a; and b; can now be found by solving 
equations, these coefficients are as below: 
Ao = —4,a, = 0, az = —2, bo = 4, b; = —8, bz = 8. 
Then, 

x(t) = —4— 227 


41 
X2(t) = 4— 8t + 8t? ay 
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Now, taking the Padé approximant of Eq. (41) ,with | = 0, m= 2, we get the solutions of the system (35): 


—4 
*110,2] (t) = 1-0.5t2 
4 
*2 [0,2 ©) ~ 442¢42¢2 
4. Results and discussion 


This section, discusses the numerical computations for predator- prey models, which have been obtained by 
the application of SAGPM. All calculations are employed by Maple 2016 software. In Tables1, 2 and 3, a 
comparison was made between the errors, the central processing unit, and the rate of convergence obtained 
using AGM, modified decomposition method (ADM) [38], He’s variational iteration method (VIM) [15], 
Homotopy perturbation method (HPM) [37], Adomian decomposition method(ADM) [39] and SAGPM. It 
was noted through the results that the proposed method is characterized by high accuracy, fewer errors, and a 
lower central processing unit when compared to other methods. Fig.1 shows the comparison between the 
solution offered by the new method SAGPM, AGM, ADM{[38], VIM [15] and exact solutions 
for x,(t) and x(t) for Example1. 


Fig.2 shows the comparison between the solution offered by the new method SAGPM, AGM, HPM [37] and 
ADM [39] for x(t) and x2(t) respectively for Example2. 


Fig.3 shows the comparison between the solution offered by the new method SAGPM, AGM, ADM [39] 
and exact solutions for x(t) and x2(t) by using second terms for Example 3. According to the calculations 
demonstrated in the tables and figures, the SAGPM processes are particularly effective at solving prey- 
predator systems, and the SAGPM is the most efficient one. It also gives high-precision solutions since it 
yields good results with solution iterations and errors. 


Table 1. Comparison of error, CPU time and rate of convergence between SAGPM, AGM,VIM, and ADM, 
when t € [0,1] for Example 1 


Functions Errors SAGPM AGM VIM [15] ADM [38] 
L2 0.06391072710 | 0.4716321433 | 0.3006935003 | 0.4716321433 
x1(t)=x2(t) Lo 0.3065 15496 1.94348450 1.343979181 1.943484504 
CPU(s) 0.015 0.031 0.031 0.015 
Rate 8.24 0.95 0.95 0.95 


Table 2. Comparison of error, CPU time and rate of convergence between SAGPM, AGM,HPM and ADM, 
when t € [0,5] for Example 2. 


Functions | Errors SAGPM AGM HPM [37] ADM [39] 

L2 0.01099356035 12.33180023 12.33178100 12.33086147 

x,(t) Lo 0.00962745319 13.41102090 13.41100000 13.41000000 
CPU(s) 0.016 0.032 0.031 0.031 
Rate 2.78 1.01 1.01 1.01 

L2 0.01797838041 | 0.1297084105 0.1296993297 | 0.1297223179 

x2(t) | 0.01226168904 | 0.1410598755 0.1410500000 | 0.1410750000 
CPU(s) 0.015 0.031 0.016 0.016 
Rate 1.22 0.55 0.55 0.55 
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Table 3. Comparison of error, CPU time and rate of convergence between SAGPM,AGM and ADM ,when 
t € [0,1] for Example 3. 


Functions Errors SAGPM AGM ADM [39] 
x1(t) L2 0.02399586585_ | 0.06355437342 | 0.06355437342 
| Fae 0.596737128 1.403262872 1.403262872 
CPU(s) 6.000 6.796 6.796 
Rate 4.175 100 =f --e--- 
xX2(t) L2 0.02226227582 | 0.2027234285 0.2027234285 
Lo 0.2586588672 3.458658867 3.458658867 
CPU(s) 6.421 7.234 7.234 
Rate 2.321 0.68 = | weeee+ 


Where; the measurement errors are defined as the following: 
Ells, = Reale = x APPTOx |2 : 
Elle. Mdina ("=e |) 


and the convergence rate is defined as the following: 


E; 
log ( Ey ) 
E 


log a) 


Where E£;_,,£; and E;,, have represented the error in steps i — 1, i, and i + 1, respectively. 


Rate = 


SAGPM 
——_ AGM 


Exact 


Fig 1. Comparison of approximate solutions obtained by SAGPM, VIM,ADM, 
and AGM with exact solution for Example 1, ( x,(t) = x2(t)). 
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| [«] . 
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03 
xl x2 
02 
01 
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t 
=== SAGPM MH HPM@® ADM=——"AGM wee SAGPM ME = =HPM Ml =ADM AGM 


Fig. 2. Comparison of approximate solutions obtained by AGM,HPM,ADM and SAGPM for x, (t) in A and 
X2(t) in B for Example 2. 


0 02 04 06 08 1 


t 
= @ AGM = ADM mm Exact = SAGPM = AGM === ADM = Exact === SAGPM 


Fig. 3. Comparison of approximate solutions obtained by AGM,ADM and SAGPM with the exact solution 
for x;(t) in A and x(t) in B for Example 3. 


5. Convergence analysis of SAGPM 


This part, illustrates how the approximate analytical solutions from SAGPM for systems (19), (27) and (35) 
converge. 
Consider the system of equations in the following form: 


xX, = F(xX4,X2), (42) 


X2 = G(X,X2), 
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where F and G are non-linear operators. The solution by using the present approach is equivalent to the 
following sequence: 


Theorem (5.1): ( Convergence of Systems )[40] 


Let H be a Hilbert space, and let F, G be an operator from H into H, and x, x2 be the exact solution of 


Eq.(42). The approximate solutions 
00 00 ti 
Dijeo™! = Lijao DI 


co er) t 
%5;,= ae 
Di si Dis a 


are convergence to exact solutionsx, , X2 respectively when 0 < p <1, 


and 


[1 )41 |] <0 []x1,|| vi ev U (0} and 30 <0 <1,|x2,,,|] so |[x2, |] ,wi en U {0}. 


Definition (5.2): For every n € N U {0},we define 


xj leaPeal| 
pp =} Jay’ ei =O ang oj =4 fei’ lei] =9 
0, otherwise 0 , otherwise 


Corollary(5.3): From the theorem above Yij2.9 X17 = Ljto G a and Yi9 X2j = Ly20 Bj a convergence 
to exact solutions x; and x2 when 0 < pj < 1and0 <g; < 1,j = 0,1,2,. 


Now, to demonstrate how the three cases' analytical approximations converge, we have applied the corollary 
as follows; 


In the first example, at t € [0,1] for x, = x2, we get 


= 4 = 0.9999875001 <1, 
X19 


2= EH = 0.5000250018 <1, 


= y= = 0.3535666492e —4 <1. 
X12 


In the second example, at t € [0,5] for x,, we get 


= 4s = 0.9939414480 <1, 
X19 


2= EH = 0.2952666150e —1 <1, 


03 = y= = 0.5302789404e —2 <1. 


and for x2 
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o, = all — 9.9939414480 <1, 


Ix2oll 


2 = = = 0.3958688457e —2 <1, 
X24 


o, = al — 9 6461243933e —3 <1. 


Ilx2all 


In the third example, at t € [0,1] for x,, we get 


=I = 0.2500052087 <1, 


P2 = t= = 0.1414184098e —3 <1, 
14 


3 = sll — 9 4166666666e —4 <1. 
[x12ll 


and for x2 


o, = eal — 9.9950455635 <1, 


Ilx2l 
2= t= = 0.7035684045e —2 <1, 
21 
r = I= = 0.3333333338e —2 <1. 
22 
6. Stability analysis 


The future state of a dynamical system can be described by its evolution, which is a fixed phenomenon. 
Stability is often the first and most crucial problem to be addressed when trying to determine how a system 
behaves. This section's goal is to look into the stability analysis of the SAGPM, which is used to solve the 
prey -predator systems. The analysis of the model is largely influenced by the equilibrium point of a system 
of nonlinear differential equations for the majority of dynamical systems. The eigenvalues of the Jacobian 
matrix of the associated dynamical system are calculated to categorize the equilibrium locations. 


In Fig.4 (A) when both x, > 0 and x2 > 0 then some trajectories are almost elliptical around the point (1,1). 
we begin with both the prey and predator population as very small. We have noticed from the figure that the 
prey first increases because there is little predation. After a while, the predator population increase because 
of a bigger food supply. This decreases the prey population as we have seen in the figure. With less 
abundance of prey, the system goes back to its initial form as the predator population starts to decline once 
more. Then the trajectory begins to repeat itself. In Fig.4 (B) We noted that the critical point is unstable. The 
solutions are unstable near point (0,0). If the critical point is stable, this leads to the attraction of non-zero 
groups to it, which leads to the extinction of both species. Therefore, the instability of this point is very 
important. However, the extinction of both species is challenging to represent because the fixed point at the 
origin is a saddle point, which is unstable. This is only possible if the whole population of prey is purposely 
wiped off, starving the predators to death. In this straightforward model, if the predators disappeared, the 
number of prey would increase unchecked. Fig.4 (C) shows some solutions approach the critical point (0,0) 
as t> 0. The critical point is stable and we naturally also see that the trajectory which corresponds to our 
initial condition also approaches the critical point. That is, the solutions are stable near this point. While the 
solutions are unstable near point (-2,3.5) because it is unstable. 


7. Conclusion 


In order to solve the predator-prey systems, a new analytical approximation method was created and used in 
this study. In comparison to the existing approaches like ADM, VIM, HPM, and AGM, the numerical results 
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from the three cases show that the new technique offers good accuracy, especially for large values of time t. 
We have observed from the results that the new method is a mathematical tool for solving predator-prey 
models with constant or variable coefficients. Also, it is a crucial technique for solving beginning and 
boundary value problems that come up in a variety of applied sciences disciplines, such as mechanics, 
physics, applied mathematics, and others. It has been noted through tables and graphs that the proposed 
method noticeable convergence and high accuracy compared to other methods. In addition, the new 
technique is an improvement of Akbari- Ganji’s method. We will focus our future work on researching a few 
fractional- order predator-prey systems. 


Pm AAAAAAAAQAAAA ! Wuneeeeeeeeeeeee 104 SSSSSSSSSAANA S| 
fe NY B NVA SSS eeerrce |} C LwwAAANNNNNNN4 VJ 
mem SSSSSAAAAAAAN 044 NVA SS ferro | SSSSSSSSSSANSN) | / 
Lf fom SSAA — NAVD Seer erorene FSSA SSSSANA § | / 
bf Zee SS SAAAASNAAAN NVbESZ Fm SSSSSSANS 4 | / 
Vs ANAAAANAAARANR \hi/ 27 SSS SSS SAAS) | 
bf Jee SS NANAAAAAANS 1 SSS SSSSSSS \ | / 
PJ ZoemS SSN RANAAAANASN i SSS | / 
PI SemeSSNANAAANAAAAAN i Dae SSSA | / 
M/ NNNANA } ————— = \/ 
1/ AANANAN ~~ | eee eee eer // \ 
ti VNAAASN JmenX 4) ff FECEECE CECE C ERT ET S| 
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Fig.4: Direction fields and trajectory to systems (19), (27), and (35) at t = 20 and using the parameter 
values a, = a2 = 0.5 and ff, = B, = 0.3 in (A), b = 0.8,c = 0.2,e = 0.5,r = 0.9 in (B), a, = —4,a2 = 
2, By = 7, Bz = 2 in (©). 
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